Report generated on 2020-02-25, 12:44 based on data in:
/mnt/research/ShadeLab/GLBRC/mapping/metaT/fullAssembly/trimStats
#HIDDEN
from IPython.display import Image
from IPython.core.display import HTML
from ipywidgets import Text, Dropdown, Output, interact
from matplotlib import pyplot as plt
from pandas import DataFrame, merge, read_csv, Series, to_datetime
import cufflinks as cf
import qgrid
cf.go_offline()
# cf.set_config_file(offline=False, world_readable=True)
cf.set_config_file(offline=True, world_readable=True)
### nbi:hide_in
GLBRC collaboration
Adina, Ashley, Nejc, Shane¶
- MetaG Sequencing and alignment Summaries
- Inputs:
- 136 Phred-33 fastq files from samples obtained from washes of Switchgrass and Miscanthus transfered from JGI. Metadata
- Illumina TruSeq version 3.2 sequencing adapters.
- Outputs:
- Inputs:
- MetaG Analysis
- Inputs:
- The MAGs assembly created by concatenating bins with $>=$ 50% commpleteness
- Functional Gene Count table, mapped to KEGG/functions using this annotation methodology. My script that follows these steps is here.
- Research Questions:
- How do functional annotations accumulate over time?
- Is time a factor that affects the functional richness of phylosphere communities (Alpha Diversity)?
- How do the samples' functional counts compare each other (PCoA)?
- Do we see changes over time in functional richness?
- Is there a difference between treatment (Fertilized vs. Unfertilized)?
- Is there a difference between host plant?
- Is there similarity between the WGS communities and the 16S data?
- How variable are metagenomes across replicate time points?
- Inputs:
- MetaT Sequencing and alignment Summaries
- Inputs:
- 74 Switchgrass samples and their respective fastq files. Metadata
- Illumina TruSeq version 3.2 sequencing adapters.
- Outputs:
- Inputs:
- MetaT Analysis
MetaG Sequencing and Alignment Summaries
Adapter Trimming and QC (Trimmomatic) Report
What is the abundance and quality of the reads in each sample?
Step 1. Split JGI files in PE files
split-paired-reads.py --gzip -1 \$sample.fastq.pe1.gz -2 \$sample.fastq.pe2.gz \$jgi_combined_fastq
Step 2. Use trimmomatic to QC reads
trimmomatic PE -phred33 \$sample.fastq.pe1.gz \$sample.fastq.pe2.gz \$sample.fastq.pe1.gz \$sample.fastq.se1.gz \$sample.fastq.pe2.gz \$sample.fastq.se2.gz \
ILLUMINACLIP:TruSeq3-PE-2.fa:2:30:10:8:TRUE \
LEADING:3 \
TRAILING:3 \
SLIDINGWINDOW:4:15 \
MINLEN:36
ILLUMINACLIP:TruSeq3-PE-2.fa:2:30:10:8:TRUE \
LEADING:3 \
TRAILING:3 \
SLIDINGWINDOW:4:15 \
MINLEN:36
Go back to the Outline
- Adapter Trimming and QC (Trimmomatic) Report
- Host Plant Alignment Report
- Fungal Alignments
- Alignment Conclusions
#HIDDEN
# <!-- # #HIDDEN
# # tImages = {"MetaG":Image("images/TrimReport_MetaG.png"),"MetaT":Image("images/MetaT_trimmomatic_plot.png")}
# # @interact
# # def trimReport(Source=["MetaG","MetaT"]):return tImages[Source]
# <img src="images/TrimReport_MetaG.png"></img>
# -->
HTML("html/MetaG_TrimStats.html")
Host Plant Alignment Report
When the reads are aligned to plant reference genomes, how many of the reads are plant-aligned reads?
Step 1. Align reads to respective host reference assembly
if [[ \$sample == *"G5"* ]]; then
SAMFILE=\$sample.SWGRASS.sam
BAMFILE=\$sample.SWGRASS.bam
bowtie2 --threads \$THREADS -x \$SWITCHGRASS
-1 \$sample.fastq.pe1.gz -2 \$sample.fastq.pe2.gz
-U \$sample.fastq.se12.gz -S \$sample.SWGRASS.sam 2>\$sample.stat
else
SAMFILE=\$sample.MISCANTS.sam
BAMFILE=\$sample.MISCANTS.bam
bowtie2 --threads \$THREADS -x \$MISCANTHUS
-1 \$sample.fastq.pe1.gz -2 \$sample.fastq.pe2.gz
-U \$sample.fastq.se12.gz -S \$sample.MISCANTS.sam 2>\$sample.stat
fi
SAMFILE=\$sample.SWGRASS.sam
BAMFILE=\$sample.SWGRASS.bam
bowtie2 --threads \$THREADS -x \$SWITCHGRASS
-1 \$sample.fastq.pe1.gz -2 \$sample.fastq.pe2.gz
-U \$sample.fastq.se12.gz -S \$sample.SWGRASS.sam 2>\$sample.stat
else
SAMFILE=\$sample.MISCANTS.sam
BAMFILE=\$sample.MISCANTS.bam
bowtie2 --threads \$THREADS -x \$MISCANTHUS
-1 \$sample.fastq.pe1.gz -2 \$sample.fastq.pe2.gz
-U \$sample.fastq.se12.gz -S \$sample.MISCANTS.sam 2>\$sample.stat
fi
Go back to the Outline
- Adapter Trimming and QC (Trimmomatic) Report
- Host Plant Alignment Report
- Fungal Alignments
- Alignment Conclusions
#HIDDEN
HTML(filename="html/BothPerc.html")
Sequencing and alignment Summaries
Fungal Alignment Report
What reads are fungal reads?
Step 1. Extract reads that don't align the plant assembly
# A. Convert sam > bam
samtools view -bS \$SAMFILE > \$BAMFILE
# B. Extract all unmapped reads (-f 4) that don't have a mate mapped (-F 256 i.e. both unmapped)
samtools view -b -f 4 -F 256 \$BAMFILE > \$sample.unmapped.bam
# C. Sort the reads
samtools sort -n \$sample.unmapped.bam -o \$sample.unmapped_sorted.bam
# D. Split the unmapped reads in R1 & R2 (back into paired end)
bedtools bamtofastq -i \$sample.unmapped_sorted.bam -fq \$sample.R1.fastq -fq2 \$sample.R2.fastq
samtools view -bS \$SAMFILE > \$BAMFILE
# B. Extract all unmapped reads (-f 4) that don't have a mate mapped (-F 256 i.e. both unmapped)
samtools view -b -f 4 -F 256 \$BAMFILE > \$sample.unmapped.bam
# C. Sort the reads
samtools sort -n \$sample.unmapped.bam -o \$sample.unmapped_sorted.bam
# D. Split the unmapped reads in R1 & R2 (back into paired end)
bedtools bamtofastq -i \$sample.unmapped_sorted.bam -fq \$sample.R1.fastq -fq2 \$sample.R2.fastq
Step 2. Align the reads to the combined fungal assemblies
bowtie2 -x \$FUNGAL -1 \$sample.R1.fastq -2 \$sample.R2.fastq -S \$sample.fungal.sam 2>\$sample.fungal.stat
Go back to the Outline
#HIDDEN
HTML(filename='html/FungalAlign.html')
Sequencing and alignment Conclusions
Sequencing Observations
What patterns are there over the season?
- The percentage of plant reads sequenced is high thoughout the season and tapers off towards the end of the season. This may be do to senescence of plant cells.
- Alignment to the 6 selected fungal assemblies is low, but picks up during the warmer months. Overall, reads aligning to the combined fungal assembly were <10% of the non-plant reads.
Annotation Summary¶
The annotations were performed using KEGG's prokaryotic peptides and eukaryotic peptides
#HIDDEN
HTML(filename="html/MAG_Stats.html")
MetaT Analysis
Overarching analysis - changes in metagenome content over time
All bins with at least 50% completeness were used in a combined assembly. The analysis below comes from the counts for each sample. These counts were normalized by the number of reads that did not align to the fungal or host assemblies.
Go back to the Outline
- Functional annotation accumulation over time
- Is time a factor that affects the functional richness of phylosphere communities (Alpha Diversity)?
- Do we see changes over time in functional richness?
- Is there a difference between treatment (Fertilized vs. Unfertilized)?
- Is there a difference between host plant?
- PCoA of functional counts
- Is there similarity between the WGS communities and the 16S data?
- How variable are metagenomes across replicate time points?
COG accumulation over time
Is time a factor that affects the functional richness of phylosphere communities (Alpha Diversity)?
|
|
|
||||||||||||||||||||||||||||||||||||
|
|
|
||||||||||||||||||||||||||||||||||||
Call: adonis(formula = dist.otu ~ map_16S\$time_numeric)
Permutation: free Number of permutations: 999
Terms added sequentially (first to last)
| Df | SumsOfSqs | MeanSqs | F.Model | R2 | Pr(>F) | |
|---|---|---|---|---|---|---|
| map_16S$time_numeric | 1 | 1.9356 | 1.93562 | 47.814 | 0.26298 | 0.001 *** |
| Residuals | 134 | 5.4246 | 0.04048 | 0.73702 | ||
| Total | 135 | 7.3603 | 1.000000 |
Call: adonis(formula = dist.otu ~ map_16S\$treatment)
Permutation: free Number of permutations: 999
Terms added sequentially (first to last)
| Df | SumsOfSqs | MeanSqs | F.Model | R2 | Pr(>F) | |
|---|---|---|---|---|---|---|
| map_16S$treatment | 1 | 0.0289 | 0.028930 | 0.52878 | 0.00393 | 0.755 |
| Residuals | 134 | 7.3313 | 0.054711 | 0.99607 | ||
| Total | 135 | 7.3603 | 1.000000 |
Call: adonis(formula = dist.otu ~ map_16S\$plant)
Permutation: free Number of permutations: 999
Terms added sequentially (first to last)
| Df | SumsOfSqs | MeanSqs | F.Model | R2 | Pr(>F) | |
|---|---|---|---|---|---|---|
| map_16S$plant | 1 | 0.1617 | 0.161704 | 3.0101 | 0.02197 | 0.027 * |
| Residuals | 134 | 7.1986 | 0.053721 | 0.97803 | ||
| Total | 135 | 7.3603 | 1.000000 |
PCoA
Go back to the Outline
- Functional annotation accumulation over time
- Is time a factor that affects the functional richness of phylosphere communities (Alpha Diversity)?
- Do we see changes over time in functional richness?
- Is there a difference between treatment (Fertilized vs. Unfertilized)?
- Is there a difference between host plant?
- PCoA of functional counts
- Is there similarity between the WGS communities and the 16S data?
- How variable are metagenomes across replicate time points?
Go back to the Outline
- Functional annotation accumulation over time
- Is time a factor that affects the functional richness of phylosphere communities (Alpha Diversity)?
- Do we see changes over time in functional richness?
- Is there a difference between treatment (Fertilized vs. Unfertilized)?
- Is there a difference between host plant?
- PCoA of functional counts
- Is there similarity between the WGS communities and the 16S data?
- How variable are metagenomes across replicate time points?
MetaT Sequencing and Alignment Summaries
What is the abundance and quality of the reads in each sample?
#HIDDEN
HTML("html/MetaT_TrimStats.html")
v1.7
Loading report..
Highlight Samples
Regex mode off
Rename Samples
Paste two columns of a tab-delimited table here (eg. from Excel).
First column should be the old name, second column the new name.
Regex mode off
Show / Hide Samples
Regex mode off
Export Plots
Save Settings
You can save the toolbox settings for this report to the browser.
Load Settings
Choose a saved report profile from the dropdown box below:
About MultiQC
This report was generated using MultiQC, version 1.7
You can see a YouTube video describing how to use MultiQC reports here: https://youtu.be/qPbIlO_KWN0
For more information about MultiQC, including other videos and extensive documentation, please visit http://multiqc.info
You can report bugs, suggest improvements and find the source code for MultiQC on GitHub: https://github.com/ewels/MultiQC
MultiQC is published in Bioinformatics:
MultiQC: Summarize analysis results for multiple tools and samples in a single report
Philip Ewels, Måns Magnusson, Sverker Lundin and Max Käller
Bioinformatics (2016)
doi: 10.1093/bioinformatics/btw354
PMID: 27312411
MultiQC v1.7
- Written by Phil Ewels,
available on GitHub.
This report uses HighCharts, jQuery, jQuery UI, Bootstrap, FileSaver.js and clipboard.js.
#HIDDEN
HTML("html/MetaT_HostAlignment.html")
v1.7
Loading report..
Highlight Samples
Regex mode off
Rename Samples
Paste two columns of a tab-delimited table here (eg. from Excel).
First column should be the old name, second column the new name.
Regex mode off
Show / Hide Samples
Regex mode off
Export Plots
Save Settings
You can save the toolbox settings for this report to the browser.
Load Settings
Choose a saved report profile from the dropdown box below:
About MultiQC
This report was generated using MultiQC, version 1.7
You can see a YouTube video describing how to use MultiQC reports here: https://youtu.be/qPbIlO_KWN0
For more information about MultiQC, including other videos and extensive documentation, please visit http://multiqc.info
You can report bugs, suggest improvements and find the source code for MultiQC on GitHub: https://github.com/ewels/MultiQC
MultiQC is published in Bioinformatics:
MultiQC: Summarize analysis results for multiple tools and samples in a single report
Philip Ewels, Måns Magnusson, Sverker Lundin and Max Käller
Bioinformatics (2016)
doi: 10.1093/bioinformatics/btw354
PMID: 27312411
Report regenerated on 2020-02-18, 21:40 based on data in:
/mnt/research/ShadeLab/GLBRC/mapping/metaT/hostRemovalFlagstats
There are 6 possible types of alignment: * PE mapped uniquely: Pair has only one occurence in the reference genome. * PE mapped discordantly uniquely: Pair has only one occurence but not in proper pair. * PE one mate mapped uniquely: One read of a pair has one occurence. * PE multimapped: Pair has multiple occurence. * PE one mate multimapped: One read of a pair has multiple occurence. * PE neither mate aligned: Pair has no occurence.
MultiQC v1.7
- Written by Phil Ewels,
available on GitHub.
This report uses HighCharts, jQuery, jQuery UI, Bootstrap, FileSaver.js and clipboard.js.
MetaT Analysis
MetaT Kallisto¶
Link to the interactive tables. Static tables shown below.
MetaT Trim Report
Link to the interactive metaT trim report.
Alignment of the MetaT reads to the MAG assembly
MetaT alignment by count
MetaT alignment by percentage
Next Steps¶
- Fix the metadata to ensure that each sample is accounted for and we have the seq files (Nejc) ✓ Completed Wednesday, Feb 19,2020. Delivered by Kiera
- Modify the R analysis script to read results from the MetaT data (Nejc & Shane)
In order for this to progress, 1 had to be complete. The next step is figure out why the metadata for the MetaT is not meshing with the R script.
- Explore the Kallisto output to determine how to extract annotations (Shane)
Kallisto has an arguement to export a bam file when it runs, but I do not know how these bam files work or if they are the same thing as a bam from bowtie2. I am working with kallisto to get a bam file now and then I will see if that file works to extract annotations from the contigs